Calculate Effect size for visits

# use ps for stool samples only
ps_full <- readRDS("/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/rds_files/ps_IMProveCF_patientsAndControls_Run1-23_18102023.rds")
# subset per material
ps_full_stool <- subset_samples(ps_full, material== "Stool")
# remove zero abundances from each dataset
ps_stool <- tax_filter(ps_full_stool, min_prevalence = 0.3) # has to be prevalent in 30% of samples

## Proportional min_prevalence given: 0.3 --> min 76/251 samples.

# vs
# calculate relative abundances
ps_stool_relab <- transform_sample_counts(ps_stool, function(x) x/sum(x))
#ps_stool_relab <-prune_taxa(taxa_sums(ps_stool_relab) > 0.45, ps_stool_relab) # > 40% prevalence over all samples, 107 most abundant ASVs

ps_stool_relab_family <- tax_glom(ps_stool_relab, "Family")
ps_stool_relab_order <- tax_glom(ps_stool_relab, "Order")
ps_stool_relab_class <- tax_glom(ps_stool_relab, "Class")
ps_stool_relab_phylum<- tax_glom(ps_stool_relab, "Phylum")

ps_stool_relab <- tax_glom(ps_stool_relab, "Genus")

ASV/Genus level

psControlV2 <-
 phyloseq::subset_samples(ps_stool_relab , visit_cal_9  %in% c("Control", "2"))

metadata_ControlV2 <- as(sample_data(psControlV2),"data.frame")

metada <- metadata_ControlV2%>%
  mutate(visit_ControlV2 = (as.numeric(visit_cal_9))-1) # Control is assigned 1, CF timepoint 0

metada<- cbind(visit_ControlV2 =metada$visit_ControlV2, subset(metada,select = -c(visit_ControlV2))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_ControlV2)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControlV2))
taxtable <- as.data.frame(tax_table(psControlV2))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  dplyr::mutate(ASV_Genus = paste(row.names(.), Genus, sep = "_")) %>%
  dplyr::select(ASV_Genus = ASV_Genus, -Genus)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9) # if I compare Control vs CF, there are no repeated measures in this model , randomVar = list("+ (1|id)", c("id"))

raw_p_ControlV2 <- metad[1]
corr_p_ControlV2 <- metad[2]
effect_size_ControlV2 <- metad[3]
status_ControlV2 <- metad[4]

psControl_V3V5 <-
 phyloseq::subset_samples(ps_stool_relab , visit_cal_cor  %in% c(NA, "3", "4", "5"))

metadata_Control_V3V5 <- as(sample_data(psControl_V3V5),"data.frame")

metadata_Control_V3V5$visit_sum

##   [1] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [10] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [19] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [28] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [37] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [46] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [55] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [64] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [73] Control Control Control Control Control Control 3-5     3-5     3-5    
##  [82] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [91] Control Control Control Control Control Control Control 3-5     3-5    
## [100] 3-5     3-5     Control Control Control Control Control Control Control
## [109] Control Control Control Control Control Control Control Control Control
## [118] Control Control Control Control
## Levels: 3-5 Control

metada <- metadata_Control_V3V5%>%
  mutate(visit_Control_V3V5 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V3V5 =metada$visit_Control_V3V5, subset(metada, select = -c(visit_Control_V3V5))) # brings visit into the first column

metada <- metada%>%
 dplyr::select(visit_Control_V3V5) # COntrols are assigned 1

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V3V5))
taxtable <- as.data.frame(tax_table(psControl_V3V5))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Genus = paste(row.names(.), Genus, sep = "_")) %>%
  dplyr::select(ASV_Genus = ASV_Genus, -Genus)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id"))

raw_pControl_V3V5 <- metad[1]
corr_pControl_V3V5 <- metad[2]
effect_sizeControl_V3V5 <- metad[3]
statusControl_V3V5 <- metad[4]

psControl_V6V7 <-
 phyloseq::subset_samples(ps_stool_relab , visit_cal_cor  %in% c(NA, "6", "7"))

metadata_Control_V6V7 <- as(sample_data(psControl_V6V7),"data.frame")

metadata_Control_V6V7$visit_cal_cor

##  [1] 6    6    7    6    7    6    6    6    6    6    6    6    6    6    6   
## [16] 6    7    6    7    7    7    6    7    7    7    6    7    7    7    6   
## [31] 7    7    7    7    7    6    7    7    <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> 7    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> 6    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [76] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 6 7

metadata_Control_V6V7$visit_sum

##  [1] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [10] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [19] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [28] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [37] 6-7     6-7     Control Control Control Control Control Control Control
## [46] Control Control Control Control Control 6-7     Control Control Control
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control 6-7     Control Control Control Control Control Control Control
## [73] Control Control Control Control Control Control Control Control Control
## [82] Control Control Control Control
## Levels: 6-7 Control

metada <- metadata_Control_V6V7%>%
  mutate(visit_Control_V6V7 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V6V7 =metada$visit_Control_V6V7, subset(metada, select = -c(visit_Control_V6V7))) # brings visit into the first column

metada <- metada%>%
dplyr::select(visit_Control_V6V7)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V6V7))
taxtable <- as.data.frame(tax_table(psControl_V6V7))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Genus = paste(row.names(.), Genus, sep = "_")) %>%
  dplyr::select(ASV_Genus = ASV_Genus, -Genus)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada, nnodes=9)#, randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V6V7 <- metad[1]
corr_p_Control_V6V7 <- metad[2]
effect_size_Control_V6V7 <- metad[3]
status_Control_V6V7 <- metad[4]

psControl_V8V9 <-
 phyloseq::subset_samples(ps_stool_relab , visit_cal_9  %in% c("Control", "8", "9"))

metadata_Control_V8V9 <- as(sample_data(psControl_V8V9),"data.frame")

metada <- metadata_Control_V8V9%>%
  mutate(visitControl_V8V9 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V8V9 =metada$visitControl_V8V9, subset(metada, select = -c(visitControl_V8V9))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V8V9)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V8V9))
taxtable <- as.data.frame(tax_table(psControl_V8V9))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Genus = paste(row.names(.), Genus, sep = "_")) %>%
  dplyr::select(ASV_Genus = ASV_Genus, -Genus)
taxtable$ASV <- row.names(taxtable)

metada$visitControl_V8V9

##  [1] 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [39] 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada, nnodes=9) #, randomVar = list("+ (1|id)", c("id")),  

raw_p_Control_V8V9 <- metad[1]
corr_p_Control_V8V9 <- metad[2]
effect_size_Control_V8V9 <- metad[3]
status_Control_V8V9 <- metad[4]

psControl_V1 <-
 phyloseq::subset_samples(ps_stool_relab , visit_cal_9  %in% c("Control", "1"))

metadata_Control_V1 <- as(sample_data(psControl_V1),"data.frame")

metada <- metadata_Control_V1%>%
  mutate(visitControl_V1 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V1 =metada$visitControl_V1, subset(metada,select = -c(visitControl_V1))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V1)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V1))
taxtable <- as.data.frame(tax_table(psControl_V1))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  

taxtable <- taxtable %>%
  mutate(ASV_Genus = paste(row.names(.), Genus, sep = "_")) %>%
  dplyr::select(ASV_Genus = ASV_Genus, -Genus)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)

raw_p_Control_V1 <- metad[1]
corr_p_Control_V1 <- metad[2]
effect_size_Control_V1 <- metad[3]
status_Control_V1 <- metad[4]

raw_p <- purrr::map2(raw_p_ControlV2,raw_pControl_V3V5,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V6V7,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V8V9,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V1,cbind)

raw_p_df <- bind_rows(raw_p)
raw_p_df <- data.frame(raw_p_df$Ps)
raw_p_df  <- raw_p_df %>%
  rownames_to_column()%>%
  mutate(p_Control_V2=visit_ControlV2)%>%
  mutate(p_Control_V3V5=visit_Control_V3V5)%>%
  mutate(p_Control_V6V7=visit_Control_V6V7)%>%
  mutate(p_Control_V8V9=visitControl_V8V9)%>%
  mutate(p_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  ##dplyr::select(-c(1:10))

corr_p <- purrr::map2(corr_p_ControlV2,corr_pControl_V3V5,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V6V7,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V8V9,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V1,cbind)

corr_p_df <- bind_rows(corr_p)
corr_p_df <- data.frame(corr_p_df$Qs)
corr_p_df  <- corr_p_df %>%
  rownames_to_column()%>%
  mutate(q_Control_V2=visit_ControlV2)%>%
  mutate(q_Control_V3V5=visit_Control_V3V5)%>%
  mutate(q_Control_V6V7=visit_Control_V6V7)%>%
  mutate(q_Control_V8V9=visitControl_V8V9)%>%
  mutate(q_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  ##dplyr::select(-c(1:10))

effect_size <- purrr::map2(effect_size_ControlV2,effect_sizeControl_V3V5,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V6V7,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V8V9,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V1,cbind)

effect_size_df <- bind_rows(effect_size)
effect_size_df <- data.frame(effect_size_df$Ds)
effect_size_df  <- effect_size_df %>%
  rownames_to_column()%>%
  mutate(d_Control_V2=visit_ControlV2)%>%
  mutate(d_Control_V3V5=visit_Control_V3V5)%>%
  mutate(d_Control_V6V7=visit_Control_V6V7)%>%
  mutate(d_Control_V8V9=visitControl_V8V9)%>%
  mutate(d_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
 # #dplyr::select(-c(1:10))


status <- purrr::map2(status_ControlV2,statusControl_V3V5,cbind)
status <- purrr::map2(status,status_Control_V6V7,cbind)
status <- purrr::map2(status,status_Control_V8V9,cbind)
status <- purrr::map2(status,status_Control_V1,cbind)

status_df <- bind_rows(status)
status_df <- data.frame(status_df$status)
status_df  <- status_df %>%
  rownames_to_column()%>%
  mutate(status_Control_V2=visit_ControlV2)%>%
  mutate(status_Control_V3V5=visit_Control_V3V5)%>%
  mutate(status_Control_V6V7=visit_Control_V6V7)%>%
  mutate(status_Control_V8V9=visitControl_V8V9)%>%
  mutate(status_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
 # #dplyr::select(-c(1:10))

effect_table <- raw_p_df%>%
  full_join(corr_p_df, by="ASV")%>%
  full_join(effect_size_df, by="ASV")%>%
  full_join(status_df, by="ASV")%>%
  full_join(taxtable, by="ASV")

effect_table_sig <- effect_table%>%
  filter(status_Control_V1=="OK_nc"|status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc")

#pivot long format

effect_table_sig_long <- effect_table_sig%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Control", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Control,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Control, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Control", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Control,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Control, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Control", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Control,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Control, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Control", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Control,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Control, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa=paste(ASV_Genus, "(Genus)"))%>%
  dplyr::select(-ASV_Genus)

effect_table_sig_long$comparison_p <- factor(effect_table_sig_long$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long%>%
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = as.factor (sign (effectSize)), size = abs (effectSize), color=(fdr))) +
  scale_shape_manual (values = c (25, 24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1","3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

# dplyr::select the entries which have OK_nc in status with all differences to Control

effect_table_sig_control <- effect_table%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc"| status_Control_V1=="OK_nc")

#pivot long format

effect_table_sig_long_control <- effect_table_sig_control%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa=paste(ASV_Genus, "(Genus)"))%>%
  dplyr::select(-ASV_Genus)

effect_table_sig_long_control$comparison_p <- factor(effect_table_sig_long_control$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_control%>%
  mutate(shape_sign= as.factor(sign(effectSize))) %>% 
  filter(shape_sign!=0) %>% #fiter out effect sizes that have no direction
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = shape_sign, size = abs (effectSize), color=(fdr)))+
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

Family level

psControlV2 <-
 phyloseq::subset_samples(ps_stool_relab_family , visit_cal_9  %in% c("Control", "V2"))

metadata_ControlV2 <- as(sample_data(psControlV2),"data.frame")

metada <- metadata_ControlV2%>%
  mutate(visit_ControlV2 = (as.numeric(visit.x))-1)

metada<- cbind(visit_ControlV2 =metada$visit_ControlV2, subset(metada, select = -c(visit_ControlV2))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_ControlV2)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControlV2))
taxtable <- as.data.frame(tax_table(psControlV2))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Family = paste(row.names(.), Family, sep = "_")) %>%
  dplyr::select(ASV_Family = ASV_Family, -Family)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,   nnodes=9)#randomVar = list("+ (1|id)", c("id")),

raw_p_ControlV2 <- metad[1]
corr_p_ControlV2 <- metad[2]
effect_size_ControlV2 <- metad[3]
status_ControlV2 <- metad[4]

psControl_V3V5 <-
 phyloseq::subset_samples(ps_stool_relab_family , visit_cal_9  %in% c("Control", "3", "4", "5"))

metadata_Control_V3V5 <- as(sample_data(psControl_V3V5),"data.frame")

metadata_Control_V3V5$visit_cal_cor

##   [1] 3    4    5    3    4    3    3    4    5    3    4    3    4    5    3   
##  [16] 4    3    4    5    3    3    5    4    3    4    3    4    5    3    4   
##  [31] 3    3    4    3    4    5    3    4    5    3    4    3    4    5    4   
##  [46] 5    3    4    5    3    4    5    5    3    5    4    3    5    4    5   
##  [61] 5    5    4    5    3    3    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
##  [76] <NA> <NA> <NA> 4    3    3    4    4    3    <NA> <NA> <NA> <NA> <NA> <NA>
##  [91] <NA> <NA> <NA> <NA> <NA> <NA> <NA> 5    5    5    4    <NA> <NA> <NA> <NA>
## [106] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [121] <NA>
## Levels: 3 4 5

metadata_Control_V3V5$visit_sum

##   [1] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [10] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [19] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [28] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [37] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [46] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [55] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [64] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [73] Control Control Control Control Control Control 3-5     3-5     3-5    
##  [82] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [91] Control Control Control Control Control Control Control 3-5     3-5    
## [100] 3-5     3-5     Control Control Control Control Control Control Control
## [109] Control Control Control Control Control Control Control Control Control
## [118] Control Control Control Control
## Levels: 3-5 Control

metada <- metadata_Control_V3V5%>%
  mutate(visit_Control_V3V5 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V3V5 =metada$visit_Control_V3V5, subset(metada,select = -c(visit_Control_V3V5))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V3V5)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V3V5))
taxtable <- as.data.frame(tax_table(psControl_V3V5))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Family = paste(row.names(.), Family, sep = "_")) %>%
  dplyr::select(ASV_Family = ASV_Family, -Family)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id"))

raw_pControl_V3V5 <- metad[1]
corr_pControl_V3V5 <- metad[2]
effect_sizeControl_V3V5 <- metad[3]
statusControl_V3V5 <- metad[4]

psControl_V6V7 <-
 phyloseq::subset_samples(ps_stool_relab_family , visit_cal_9  %in% c("Control", "6", "7"))

metadata_Control_V6V7 <- as(sample_data(psControl_V6V7),"data.frame")

metadata_Control_V6V7$visit_cal_cor

##  [1] 6    6    7    6    7    6    6    6    6    6    6    6    6    6    6   
## [16] 6    7    6    7    7    7    6    7    7    7    6    7    7    7    6   
## [31] 7    7    7    7    7    6    7    7    <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> 7    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> 6    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [76] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 6 7

metadata_Control_V6V7$visit_sum

##  [1] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [10] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [19] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [28] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [37] 6-7     6-7     Control Control Control Control Control Control Control
## [46] Control Control Control Control Control 6-7     Control Control Control
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control 6-7     Control Control Control Control Control Control Control
## [73] Control Control Control Control Control Control Control Control Control
## [82] Control Control Control Control
## Levels: 6-7 Control

metada <- metadata_Control_V6V7%>%
  mutate(visit_Control_V6V7 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V6V7 =metada$visit_Control_V6V7, subset(metada,select = -c(visit_Control_V6V7))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V6V7)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V6V7))
taxtable <- as.data.frame(tax_table(psControl_V6V7))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Family = paste(row.names(.), Family, sep = "_")) %>%
  dplyr::select(ASV_Family = ASV_Family, -Family)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#, randomVar = list("+ (1|id)", c("id")),

raw_p_Control_V6V7 <- metad[1]
corr_p_Control_V6V7 <- metad[2]
effect_size_Control_V6V7 <- metad[3]
status_Control_V6V7 <- metad[4]

psControl_V8V9 <-
 phyloseq::subset_samples(ps_stool_relab_family, visit_cal_cor  %in% c(NA, "8", "9"))

metadata_Control_V8V9 <- as(sample_data(psControl_V8V9),"data.frame")

metadata_Control_V8V9$visit_cal_cor

##  [1] 8    8    8    8    8    8    8    8    8    8    8    <NA> <NA> <NA> <NA>
## [16] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    9    8    9    9    8    9   
## [31] 9    8    9    9    9    9    8    8    9    9    <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 8 9

metadata_Control_V8V9$visit_sum

##  [1] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [10] 8-10    8-10    Control Control Control Control Control Control Control
## [19] Control Control Control Control Control 8-10    8-10    8-10    8-10   
## [28] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [37] 8-10    8-10    8-10    8-10    Control Control Control Control Control
## [46] Control Control Control Control Control Control Control Control 8-10   
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control Control Control Control Control Control Control Control Control
## [73] Control Control
## Levels: 8-10 Control

metada <- metadata_Control_V8V9%>%
  mutate(visitControl_V8V9 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V8V9 =metada$visitControl_V8V9, subset(metada,select = -c(visitControl_V8V9))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V8V9, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V8V9))
taxtable <- as.data.frame(tax_table(psControl_V8V9))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Family = paste(row.names(.), Family, sep = "_")) %>%
  dplyr::select(ASV_Family = ASV_Family, -Family)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,   nnodes=9)

#randomVar = list("+ (1|id)", c("id")),
raw_p_Control_V8V9 <- metad[1]
corr_p_Control_V8V9 <- metad[2]
effect_size_Control_V8V9 <- metad[3]
status_Control_V8V9 <- metad[4]

psControl_V1 <-
 phyloseq::subset_samples(ps_stool_relab_family , visit_cal_9  %in% c("Control", "1"))

metadata_Control_Control <- as(sample_data(psControl_V1),"data.frame")

metada <- metadata_Control_Control%>%
  mutate(visitControl_V1 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V1 =metada$visitControl_V1, subset(metada,select = -c(visitControl_V1))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V1)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V1))
taxtable <- as.data.frame(tax_table(psControl_V1))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable %>%
  mutate(ASV_Family = paste(row.names(.), Family, sep = "_")) %>%
  dplyr::select(ASV_Family = ASV_Family, -Family)
taxtable$ASV <- row.names(taxtable)

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)

raw_p_Control_V1 <- metad[1]
corr_p_Control_V1 <- metad[2]
effect_size_Control_V1 <- metad[3]
status_Control_V1 <- metad[4]

raw_p <- purrr::map2(raw_p_ControlV2,raw_pControl_V3V5,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V6V7,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V8V9,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V1,cbind)

raw_p_df <- bind_rows(raw_p)
raw_p_df <- data.frame(raw_p_df$Ps)
raw_p_df  <- raw_p_df %>%
  rownames_to_column()%>%
  mutate(p_Control_V2=visit_ControlV2)%>%
  mutate(p_Control_V3V5=visit_Control_V3V5)%>%
  mutate(p_Control_V6V7=visit_Control_V6V7)%>%
  mutate(p_Control_V8V9=visitControl_V8V9)%>%
  mutate(p_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

corr_p <- purrr::map2(corr_p_ControlV2,corr_pControl_V3V5,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V6V7,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V8V9,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V1,cbind)

corr_p_df <- bind_rows(corr_p)
corr_p_df <- data.frame(corr_p_df$Qs)
corr_p_df  <- corr_p_df %>%
  rownames_to_column()%>%
  mutate(q_Control_V2=visit_ControlV2)%>%
  mutate(q_Control_V3V5=visit_Control_V3V5)%>%
  mutate(q_Control_V6V7=visit_Control_V6V7)%>%
  mutate(q_Control_V8V9=visitControl_V8V9)%>%
  mutate(q_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_size <- purrr::map2(effect_size_ControlV2,effect_sizeControl_V3V5,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V6V7,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V8V9,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V1,cbind)

effect_size_df <- bind_rows(effect_size)
effect_size_df <- data.frame(effect_size_df$Ds)
effect_size_df  <- effect_size_df %>%
  rownames_to_column()%>%
  mutate(d_Control_V2=visit_ControlV2)%>%
  mutate(d_Control_V3V5=visit_Control_V3V5)%>%
  mutate(d_Control_V6V7=visit_Control_V6V7)%>%
  mutate(d_Control_V8V9=visitControl_V8V9)%>%
  mutate(d_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

status <- purrr::map2(status_ControlV2,statusControl_V3V5,cbind)
status <- purrr::map2(status,status_Control_V6V7,cbind)
status <- purrr::map2(status,status_Control_V8V9,cbind)
status <- purrr::map2(status,status_Control_V1,cbind)

status_df <- bind_rows(status)
status_df <- data.frame(status_df$status)
status_df  <- status_df %>%
  rownames_to_column()%>%
  mutate(status_Control_V2=visit_ControlV2)%>%
  mutate(status_Control_V3V5=visit_Control_V3V5)%>%
  mutate(status_Control_V6V7=visit_Control_V6V7)%>%
  mutate(status_Control_V8V9=visitControl_V8V9)%>%
  mutate(status_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_table_family <- raw_p_df%>%
  full_join(corr_p_df, by="ASV")%>%
  full_join(effect_size_df, by="ASV")%>%
  full_join(status_df, by="ASV")%>%
  full_join(taxtable, by="ASV")

# dplyr::select the entries which have OK_nc in status
effect_table_sig_family <- effect_table_family%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc")

#pivot long format

effect_table_sig_long_family <- effect_table_sig_family%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa=paste(ASV_Family, "(Family)"))%>%
  dplyr::select(-ASV_Family)

effect_table_sig_long_family$comparison_p <- factor(effect_table_sig_long_family$comparison_p, levels = c("Control_V1"," Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_family%>%
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = as.factor (sign (effectSize)), size = abs (effectSize), color=(fdr))) +
  scale_shape_manual (values = c (25, 24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

# dplyr::select the entries which have OK_nc in status with all differences to Control

effect_table_sig_control <- effect_table_family%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc"| status_Control_V1=="OK_nc")

#pivot long format

effect_table_sig_long_control_family <- effect_table_sig_control%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa=paste(ASV_Family, "(Family)"))%>%
  dplyr::select(-ASV_Family)

effect_table_sig_long_control_family$comparison_p <- factor(effect_table_sig_long_control_family$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_control_family%>%
  mutate(shape_sign= as.factor(sign(effectSize))) %>% 
  filter(shape_sign!=0) %>% #fiter out effect sizes that have no direction
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = shape_sign, size = abs (effectSize), color=(fdr)))+
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("3", "6-12", "15-18", "21-24", "Control"))+
  labs(x= "months from treatment start")

Order level

psControlV2 <-
 phyloseq::subset_samples(ps_stool_relab_order , visit_cal_cor  %in% c(NA, "2"))

metadata_ControlV2 <- as(sample_data(psControlV2),"data.frame")

metada <- metadata_ControlV2%>%
  mutate(visit_ControlV2 = (as.numeric(visit.x))-1) %>% 
  mutate(visit_ControlV2 = case_when(is.na(visit_ControlV2)~1, T~0))
metada$visit_ControlV2

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

metada<- cbind(visit_ControlV2 =metada$visit_ControlV2, subset(metada, select = -c(visit_ControlV2))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_ControlV2)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControlV2))
taxtable <- as.data.frame(tax_table(psControlV2))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Order)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#, randomVar = list("+ (1|id)", c("id")),

raw_p_ControlV2 <- metad[1]
corr_p_ControlV2 <- metad[2]
effect_size_ControlV2 <- metad[3]
status_ControlV2 <- metad[4]

psControl_V3V5 <-
 phyloseq::subset_samples(ps_stool_relab_order , visit_cal_cor  %in% c(NA, "3", "4", "5"))

metadata_Control_V3V5 <- as(sample_data(psControl_V3V5),"data.frame")

metadata_Control_V3V5$visit_cal_cor

##   [1] 3    4    5    3    4    3    3    4    5    3    4    3    4    5    3   
##  [16] 4    3    4    5    3    3    5    4    3    4    3    4    5    3    4   
##  [31] 3    3    4    3    4    5    3    4    5    3    4    3    4    5    4   
##  [46] 5    3    4    5    3    4    5    5    3    5    4    3    5    4    5   
##  [61] 5    5    4    5    3    3    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
##  [76] <NA> <NA> <NA> 4    3    3    4    4    3    <NA> <NA> <NA> <NA> <NA> <NA>
##  [91] <NA> <NA> <NA> <NA> <NA> <NA> <NA> 5    5    5    4    <NA> <NA> <NA> <NA>
## [106] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [121] <NA>
## Levels: 3 4 5

metadata_Control_V3V5$visit_sum

##   [1] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [10] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [19] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [28] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [37] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [46] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [55] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [64] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [73] Control Control Control Control Control Control 3-5     3-5     3-5    
##  [82] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [91] Control Control Control Control Control Control Control 3-5     3-5    
## [100] 3-5     3-5     Control Control Control Control Control Control Control
## [109] Control Control Control Control Control Control Control Control Control
## [118] Control Control Control Control
## Levels: 3-5 Control

metada <- metadata_Control_V3V5%>%
  mutate(visit_Control_V3V5 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V3V5 =metada$visit_Control_V3V5, subset(metada, select = -c(visit_Control_V3V5))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V3V5)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V3V5))
taxtable <- as.data.frame(tax_table(psControl_V3V5))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Order)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada, nnodes=9) #randomVar = list("+ (1|id)", c("id")), 

raw_pControl_V3V5 <- metad[1]
corr_pControl_V3V5 <- metad[2]
effect_sizeControl_V3V5 <- metad[3]
statusControl_V3V5 <- metad[4]

psControl_V6V7 <-
 phyloseq::subset_samples(ps_stool_relab_order , visit_cal_cor  %in% c(NA, "6", "7"))

metadata_Control_V6V7 <- as(sample_data(psControl_V6V7),"data.frame")

metadata_Control_V6V7$visit_cal_cor

##  [1] 6    6    7    6    7    6    6    6    6    6    6    6    6    6    6   
## [16] 6    7    6    7    7    7    6    7    7    7    6    7    7    7    6   
## [31] 7    7    7    7    7    6    7    7    <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> 7    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> 6    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [76] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 6 7

metadata_Control_V6V7$visit_sum

##  [1] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [10] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [19] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [28] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [37] 6-7     6-7     Control Control Control Control Control Control Control
## [46] Control Control Control Control Control 6-7     Control Control Control
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control 6-7     Control Control Control Control Control Control Control
## [73] Control Control Control Control Control Control Control Control Control
## [82] Control Control Control Control
## Levels: 6-7 Control

metada <- metadata_Control_V6V7%>%
  mutate(visit_Control_V6V7 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V6V7 =metada$visit_Control_V6V7, subset(metada, select = -c(visit_Control_V6V7))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V6V7)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V6V7))
taxtable <- as.data.frame(tax_table(psControl_V6V7))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Order)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9) #randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V6V7 <- metad[1]
corr_p_Control_V6V7 <- metad[2]
effect_size_Control_V6V7 <- metad[3]
status_Control_V6V7 <- metad[4]

psControl_V8V9 <-
 phyloseq::subset_samples(ps_stool_relab_order, visit_cal_9  %in% c("Control", "8", "9"))

metadata_Control_V8V9 <- as(sample_data(psControl_V8V9),"data.frame")

metadata_Control_V8V9$visit_cal_cor

##  [1] 8    8    8    8    8    8    8    8    8    8    8    <NA> <NA> <NA> <NA>
## [16] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    10   9    8    9    9    8   
## [31] 9    9    8    9    9    9    9    8    8    9    9    <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 8 9 10

metadata_Control_V8V9$visit_sum

##  [1] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [10] 8-10    8-10    Control Control Control Control Control Control Control
## [19] Control Control Control Control Control 8-10    8-10    8-10    8-10   
## [28] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [37] 8-10    8-10    8-10    8-10    8-10    Control Control Control Control
## [46] Control Control Control Control Control Control Control Control Control
## [55] 8-10    Control Control Control Control Control Control Control Control
## [64] Control Control Control Control Control Control Control Control Control
## [73] Control Control Control
## Levels: 8-10 Control

metada <- metadata_Control_V8V9%>%
  mutate(visitControl_V8V9 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V8V9 =metada$visitControl_V8V9, subset(metada, select = -c(visitControl_V8V9))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V8V9)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V8V9))
taxtable <- as.data.frame(tax_table(psControl_V8V9))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Order)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V8V9 <- metad[1]
corr_p_Control_V8V9 <- metad[2]
effect_size_Control_V8V9 <- metad[3]
status_Control_V8V9 <- metad[4]

psControl_V1 <-
 phyloseq::subset_samples(ps_stool_relab_order , visit_cal_9  %in% c("Control", "1"))

metadata_Control_Control <- as(sample_data(psControl_V1),"data.frame")

metada <- metadata_Control_Control%>%
  mutate(visitControl_V1 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V1 =metada$visitControl_V1, subset(metada, select = -c(visitControl_V1))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V1)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V1))
taxtable <- as.data.frame(tax_table(psControl_V1))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Order)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)

raw_p_Control_V1 <- metad[1]
corr_p_Control_V1 <- metad[2]
effect_size_Control_V1 <- metad[3]
status_Control_V1 <- metad[4]

raw_p <- purrr::map2(raw_p_ControlV2,raw_pControl_V3V5,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V6V7,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V8V9,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V1,cbind)

raw_p_df <- bind_rows(raw_p)
raw_p_df <- data.frame(raw_p_df$Ps)
raw_p_df  <- raw_p_df %>%
  rownames_to_column()%>%
  mutate(p_Control_V2=visit_ControlV2)%>%
  mutate(p_Control_V3V5=visit_Control_V3V5)%>%
  mutate(p_Control_V6V7=visit_Control_V6V7)%>%
  mutate(p_Control_V8V9=visitControl_V8V9)%>%
  mutate(p_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

corr_p <- purrr::map2(corr_p_ControlV2,corr_pControl_V3V5,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V6V7,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V8V9,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V1,cbind)

corr_p_df <- bind_rows(corr_p)
corr_p_df <- data.frame(corr_p_df$Qs)
corr_p_df  <- corr_p_df %>%
  rownames_to_column()%>%
  mutate(q_Control_V2=visit_ControlV2)%>%
  mutate(q_Control_V3V5=visit_Control_V3V5)%>%
  mutate(q_Control_V6V7=visit_Control_V6V7)%>%
  mutate(q_Control_V8V9=visitControl_V8V9)%>%
  mutate(q_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_size <- purrr::map2(effect_size_ControlV2,effect_sizeControl_V3V5,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V6V7,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V8V9,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V1,cbind)

effect_size_df <- bind_rows(effect_size)
effect_size_df <- data.frame(effect_size_df$Ds)
effect_size_df  <- effect_size_df %>%
  rownames_to_column()%>%
  mutate(d_Control_V2=visit_ControlV2)%>%
  mutate(d_Control_V3V5=visit_Control_V3V5)%>%
  mutate(d_Control_V6V7=visit_Control_V6V7)%>%
  mutate(d_Control_V8V9=visitControl_V8V9)%>%
  mutate(d_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

status <- purrr::map2(status_ControlV2,statusControl_V3V5,cbind)
status <- purrr::map2(status,status_Control_V6V7,cbind)
status <- purrr::map2(status,status_Control_V8V9,cbind)
status <- purrr::map2(status,status_Control_V1,cbind)

status_df <- bind_rows(status)
status_df <- data.frame(status_df$status)
status_df  <- status_df %>%
  rownames_to_column()%>%
  mutate(status_Control_V2=visit_ControlV2)%>%
  mutate(status_Control_V3V5=visit_Control_V3V5)%>%
  mutate(status_Control_V6V7=visit_Control_V6V7)%>%
  mutate(status_Control_V8V9=visitControl_V8V9)%>%
  mutate(status_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_table_order <- raw_p_df%>%
  full_join(corr_p_df, by="ASV")%>%
  full_join(effect_size_df, by="ASV")%>%
  full_join(status_df, by="ASV")%>%
  full_join(taxtable, by="ASV")

# dplyr::select the entries which have OK_nc in status
effect_table_sig_order <- effect_table_order%>%
  filter(status_Control_V1=="OK_nc"| status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc")

#pivot long format

effect_table_sig_long_order <- effect_table_sig_order%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Order, "(Order)"))%>%
  dplyr::select(-Order)

effect_table_sig_long_order$comparison_p <- factor(effect_table_sig_long_order$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_order%>%
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = as.factor (sign (effectSize)), size = abs (effectSize), color=(fdr))) +
  scale_shape_manual (values = c (25, 24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

# dplyr::select the entries which have OK_nc in status with all differences to Control

effect_table_sig_control <- effect_table_order%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc"| status_Control_V1=="OK_nc")

#pivot long format

effect_table_sig_long_control_order<- effect_table_sig_control%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Order, "(Order)"))%>%
  dplyr::select(-Order)

effect_table_sig_long_control_order$comparison_p <- factor(effect_table_sig_long_control_order$comparison_p, levels = c("Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9","Control_Control"))

effect_table_sig_long_control_order%>%
  mutate(shape_sign= as.factor(sign(effectSize))) %>% 
  filter(shape_sign!=0) %>% #fiter out effect sizes that have no direction
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = shape_sign, size = abs (effectSize), color=(fdr)))+
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("3", "6-12", "15-18", "21-24", "Control"))+
  labs(x= "months from treatment start")

class level

psControlV2 <-
 phyloseq::subset_samples(ps_stool_relab_class , visit_cal_cor  %in% c(NA, "2"))

metadata_ControlV2 <- as(sample_data(psControlV2),"data.frame")

metada <- metadata_ControlV2%>%
  mutate(visit_ControlV2 = (as.numeric(visit.x))-1) %>% 
  mutate(visit_ControlV2 = case_when(is.na(visit_ControlV2)~1, T~0))

metada<- cbind(visit_ControlV2 =metada$visit_ControlV2, subset(metada, select = -c(visit_ControlV2))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_ControlV2)

metada$visit_ControlV2

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1
## [39] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

# run metadeconfoundR
features <- as.data.frame(otu_table(psControlV2))
taxtable <- as.data.frame(tax_table(psControlV2))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Class)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_ControlV2 <- metad[1]
corr_p_ControlV2 <- metad[2]
effect_size_ControlV2 <- metad[3]
status_ControlV2 <- metad[4]

psControl_V3V5 <-
 phyloseq::subset_samples(ps_stool_relab_class , visit_cal_cor  %in% c(NA, "3", "4", "5"))

metadata_Control_V3V5 <- as(sample_data(psControl_V3V5),"data.frame")

metadata_Control_V3V5$visit_cal_cor

##   [1] 3    4    5    3    4    3    3    4    5    3    4    3    4    5    3   
##  [16] 4    3    4    5    3    3    5    4    3    4    3    4    5    3    4   
##  [31] 3    3    4    3    4    5    3    4    5    3    4    3    4    5    4   
##  [46] 5    3    4    5    3    4    5    5    3    5    4    3    5    4    5   
##  [61] 5    5    4    5    3    3    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
##  [76] <NA> <NA> <NA> 4    3    3    4    4    3    <NA> <NA> <NA> <NA> <NA> <NA>
##  [91] <NA> <NA> <NA> <NA> <NA> <NA> <NA> 5    5    5    4    <NA> <NA> <NA> <NA>
## [106] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [121] <NA>
## Levels: 3 4 5

metadata_Control_V3V5$visit_sum

##   [1] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [10] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [19] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [28] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [37] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [46] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [55] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [64] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [73] Control Control Control Control Control Control 3-5     3-5     3-5    
##  [82] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [91] Control Control Control Control Control Control Control 3-5     3-5    
## [100] 3-5     3-5     Control Control Control Control Control Control Control
## [109] Control Control Control Control Control Control Control Control Control
## [118] Control Control Control Control
## Levels: 3-5 Control

metada <- metadata_Control_V3V5%>%
  mutate(visit_Control_V3V5 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V3V5 =metada$visit_Control_V3V5, subset(metada, select = -c(visit_Control_V3V5))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V3V5)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V3V5))
taxtable <- as.data.frame(tax_table(psControl_V3V5))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Class)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_pControl_V3V5 <- metad[1]
corr_pControl_V3V5 <- metad[2]
effect_sizeControl_V3V5 <- metad[3]
statusControl_V3V5 <- metad[4]

psControl_V6V7 <-
 phyloseq::subset_samples(ps_stool_relab_class , visit_cal_cor  %in% c(NA, "6", "7"))

metadata_Control_V6V7 <- as(sample_data(psControl_V6V7),"data.frame")

metadata_Control_V6V7$visit_cal_cor

##  [1] 6    6    7    6    7    6    6    6    6    6    6    6    6    6    6   
## [16] 6    7    6    7    7    7    6    7    7    7    6    7    7    7    6   
## [31] 7    7    7    7    7    6    7    7    <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> 7    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> 6    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [76] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 6 7

metadata_Control_V6V7$visit_sum

##  [1] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [10] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [19] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [28] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [37] 6-7     6-7     Control Control Control Control Control Control Control
## [46] Control Control Control Control Control 6-7     Control Control Control
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control 6-7     Control Control Control Control Control Control Control
## [73] Control Control Control Control Control Control Control Control Control
## [82] Control Control Control Control
## Levels: 6-7 Control

metada <- metadata_Control_V6V7%>%
  mutate(visit_Control_V6V7 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V6V7 =metada$visit_Control_V6V7, subset(metada, select = -c(visit_Control_V6V7))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V6V7, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V6V7))
taxtable <- as.data.frame(tax_table(psControl_V6V7))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Class)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,   nnodes=9)#randomVar =  c("id"),

raw_p_Control_V6V7 <- metad[1]
corr_p_Control_V6V7 <- metad[2]
effect_size_Control_V6V7 <- metad[3]
status_Control_V6V7 <- metad[4]

psControl_V8V9 <-
 phyloseq::subset_samples(ps_stool_relab_class, visit_cal_cor  %in% c(NA, "8", "9","10"))

metadata_Control_V8V9 <- as(sample_data(psControl_V8V9),"data.frame")

metadata_Control_V8V9$visit_cal_cor

##  [1] 8    8    8    8    8    8    8    8    8    8    8    <NA> <NA> <NA> <NA>
## [16] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    10   9    8    9    9    8   
## [31] 9    9    8    9    9    9    9    8    8    9    9    <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 8 9 10

metadata_Control_V8V9$visit_sum

##  [1] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [10] 8-10    8-10    Control Control Control Control Control Control Control
## [19] Control Control Control Control Control 8-10    8-10    8-10    8-10   
## [28] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [37] 8-10    8-10    8-10    8-10    8-10    Control Control Control Control
## [46] Control Control Control Control Control Control Control Control Control
## [55] 8-10    Control Control Control Control Control Control Control Control
## [64] Control Control Control Control Control Control Control Control Control
## [73] Control Control Control
## Levels: 8-10 Control

metada <- metadata_Control_V8V9%>%
  mutate(visitControl_V8V9 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V8V9 =metada$visitControl_V8V9, subset(metada, select = -c(visitControl_V8V9))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V8V9)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V8V9))
taxtable <- as.data.frame(tax_table(psControl_V8V9))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Class)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V8V9 <- metad[1]
corr_p_Control_V8V9 <- metad[2]
effect_size_Control_V8V9 <- metad[3]
status_Control_V8V9 <- metad[4]

psControl_V1 <-
 phyloseq::subset_samples(ps_stool_relab_class , visit_cal_9  %in% c("Control", "1"))

metadata_Control_Control <- as(sample_data(psControl_V1),"data.frame")

metada <- metadata_Control_Control%>%
  mutate(visitControl_V1 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V1 =metada$visitControl_V1, subset(metada, select = -c(visitControl_V1))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V1)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V1))
taxtable <- as.data.frame(tax_table(psControl_V1))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Class)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)

raw_p_Control_V1 <- metad[1]
corr_p_Control_V1 <- metad[2]
effect_size_Control_V1 <- metad[3]
status_Control_V1 <- metad[4]

raw_p <- purrr::map2(raw_p_ControlV2,raw_pControl_V3V5,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V6V7,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V8V9,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V1,cbind)

raw_p_df <- bind_rows(raw_p)
raw_p_df <- data.frame(raw_p_df$Ps)
raw_p_df  <- raw_p_df %>%
  rownames_to_column()%>%
  mutate(p_Control_V2=visit_ControlV2)%>%
  mutate(p_Control_V3V5=visit_Control_V3V5)%>%
  mutate(p_Control_V6V7=visit_Control_V6V7)%>%
  mutate(p_Control_V8V9=visitControl_V8V9)%>%
  mutate(p_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

corr_p <- purrr::map2(corr_p_ControlV2,corr_pControl_V3V5,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V6V7,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V8V9,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V1,cbind)

corr_p_df <- bind_rows(corr_p)
corr_p_df <- data.frame(corr_p_df$Qs)
corr_p_df  <- corr_p_df %>%
  rownames_to_column()%>%
  mutate(q_Control_V2=visit_ControlV2)%>%
  mutate(q_Control_V3V5=visit_Control_V3V5)%>%
  mutate(q_Control_V6V7=visit_Control_V6V7)%>%
  mutate(q_Control_V8V9=visitControl_V8V9)%>%
  mutate(q_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_size <- purrr::map2(effect_size_ControlV2,effect_sizeControl_V3V5,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V6V7,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V8V9,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V1,cbind)

effect_size_df <- bind_rows(effect_size)
effect_size_df <- data.frame(effect_size_df$Ds)
effect_size_df  <- effect_size_df %>%
  rownames_to_column()%>%
  mutate(d_Control_V2=visit_ControlV2)%>%
  mutate(d_Control_V3V5=visit_Control_V3V5)%>%
  mutate(d_Control_V6V7=visit_Control_V6V7)%>%
  mutate(d_Control_V8V9=visitControl_V8V9)%>%
  mutate(d_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

status <- purrr::map2(status_ControlV2,statusControl_V3V5,cbind)
status <- purrr::map2(status,status_Control_V6V7,cbind)
status <- purrr::map2(status,status_Control_V8V9,cbind)
status <- purrr::map2(status,status_Control_V1,cbind)

status_df <- bind_rows(status)
status_df <- data.frame(status_df$status)
status_df  <- status_df %>%
  rownames_to_column()%>%
  mutate(status_Control_V2=visit_ControlV2)%>%
  mutate(status_Control_V3V5=visit_Control_V3V5)%>%
  mutate(status_Control_V6V7=visit_Control_V6V7)%>%
  mutate(status_Control_V8V9=visitControl_V8V9)%>%
  mutate(status_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_table_class <- raw_p_df%>%
  full_join(corr_p_df, by="ASV")%>%
  full_join(effect_size_df, by="ASV")%>%
  full_join(status_df, by="ASV")%>%
  full_join(taxtable, by="ASV")

# dplyr::select the entries which have OK_nc in status
effect_table_sig_class <- effect_table_class%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc")

#pivot long format

effect_table_sig_long_class <- effect_table_sig_class%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Class, "(Class)"))%>%
  dplyr::select(-Class)

effect_table_sig_long_class$comparison_p <- factor(effect_table_sig_long_class$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_class%>%
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = as.factor (sign (effectSize)), size = abs (effectSize), color=(fdr))) +
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

# dplyr::select the entries which have OK_nc in status with all differences to Control

effect_table_sig_control <- effect_table_class%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc"| status_Control_V1=="OK_nc")

#pivot long format

effect_table_sig_long_control_class<- effect_table_sig_control%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Class, "(Class)"))%>%
  dplyr::select(-Class)

effect_table_sig_long_control_class$comparison_p <- factor(effect_table_sig_long_control_class$comparison_p, levels = c("Control_V1", "Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_control_class%>%
  mutate(shape_sign= as.factor(sign(effectSize))) %>% 
  filter(shape_sign!=0) %>% #fiter out effect sizes that have no direction
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = shape_sign, size = abs (effectSize), color=(fdr)))+
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

phylum level

psControlV2 <-
 phyloseq::subset_samples(ps_stool_relab_phylum , visit_cal_cor  %in% c(NA, "2"))

metadata_ControlV2 <- as(sample_data(psControlV2),"data.frame")

metada <- metadata_ControlV2%>%
  mutate(visit_ControlV2 = (as.numeric(visit.x))-1) %>% 
  mutate(visit_ControlV2 = case_when(is.na(visit_ControlV2)~1, T~0))

metada<- cbind(visit_ControlV2 =metada$visit_ControlV2, subset(metada, select = -c(visit_ControlV2))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_ControlV2)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControlV2))
taxtable <- as.data.frame(tax_table(psControlV2))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Phylum)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_ControlV2 <- metad[1]
corr_p_ControlV2 <- metad[2]
effect_size_ControlV2 <- metad[3]
status_ControlV2 <- metad[4]

psControl_V3V5 <-
 phyloseq::subset_samples(ps_stool_relab_phylum , visit_cal_cor  %in% c(NA, "3", "4", "5"))

metadata_Control_V3V5 <- as(sample_data(psControl_V3V5),"data.frame")

metadata_Control_V3V5$visit_cal_cor

##   [1] 3    4    5    3    4    3    3    4    5    3    4    3    4    5    3   
##  [16] 4    3    4    5    3    3    5    4    3    4    3    4    5    3    4   
##  [31] 3    3    4    3    4    5    3    4    5    3    4    3    4    5    4   
##  [46] 5    3    4    5    3    4    5    5    3    5    4    3    5    4    5   
##  [61] 5    5    4    5    3    3    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
##  [76] <NA> <NA> <NA> 4    3    3    4    4    3    <NA> <NA> <NA> <NA> <NA> <NA>
##  [91] <NA> <NA> <NA> <NA> <NA> <NA> <NA> 5    5    5    4    <NA> <NA> <NA> <NA>
## [106] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [121] <NA>
## Levels: 3 4 5

metadata_Control_V3V5$visit_sum

##   [1] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [10] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [19] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [28] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [37] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [46] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [55] 3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5     3-5    
##  [64] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [73] Control Control Control Control Control Control 3-5     3-5     3-5    
##  [82] 3-5     3-5     3-5     Control Control Control Control Control Control
##  [91] Control Control Control Control Control Control Control 3-5     3-5    
## [100] 3-5     3-5     Control Control Control Control Control Control Control
## [109] Control Control Control Control Control Control Control Control Control
## [118] Control Control Control Control
## Levels: 3-5 Control

metada <- metadata_Control_V3V5%>%
  mutate(visit_Control_V3V5 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V3V5 =metada$visit_Control_V3V5, subset(metada, select = -c(visit_Control_V3V5))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V3V5, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V3V5))
taxtable <- as.data.frame(tax_table(psControl_V3V5))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Phylum)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_pControl_V3V5 <- metad[1]
corr_pControl_V3V5 <- metad[2]
effect_sizeControl_V3V5 <- metad[3]
statusControl_V3V5 <- metad[4]

psControl_V6V7 <-
 phyloseq::subset_samples(ps_stool_relab_phylum , visit_cal_cor  %in% c(NA, "6", "7"))

metadata_Control_V6V7 <- as(sample_data(psControl_V6V7),"data.frame")

metadata_Control_V6V7$visit_cal_cor

##  [1] 6    6    7    6    7    6    6    6    6    6    6    6    6    6    6   
## [16] 6    7    6    7    7    7    6    7    7    7    6    7    7    7    6   
## [31] 7    7    7    7    7    6    7    7    <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> 7    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> 6    <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## [76] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 6 7

metadata_Control_V6V7$visit_sum

##  [1] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [10] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [19] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [28] 6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7     6-7    
## [37] 6-7     6-7     Control Control Control Control Control Control Control
## [46] Control Control Control Control Control 6-7     Control Control Control
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control 6-7     Control Control Control Control Control Control Control
## [73] Control Control Control Control Control Control Control Control Control
## [82] Control Control Control Control
## Levels: 6-7 Control

metada <- metadata_Control_V6V7%>%
  mutate(visit_Control_V6V7 = (as.numeric(visit_sum))-1)

metada<- cbind(visit_Control_V6V7 =metada$visit_Control_V6V7, subset(metada,select = -c(visit_Control_V6V7))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visit_Control_V6V7, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V6V7))
taxtable <- as.data.frame(tax_table(psControl_V6V7))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Phylum)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable,select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V6V7 <- metad[1]
corr_p_Control_V6V7 <- metad[2]
effect_size_Control_V6V7 <- metad[3]
status_Control_V6V7 <- metad[4]

psControl_V8V9 <-
 phyloseq::subset_samples(ps_stool_relab_phylum, visit_cal_cor  %in% c(NA, "8", "9"))

metadata_Control_V8V9 <- as(sample_data(psControl_V8V9),"data.frame")

metadata_Control_V8V9$visit_cal_cor

##  [1] 8    8    8    8    8    8    8    8    8    8    8    <NA> <NA> <NA> <NA>
## [16] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    9    8    9    9    8    9   
## [31] 9    8    9    9    9    9    8    8    9    9    <NA> <NA> <NA> <NA> <NA>
## [46] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> 9    <NA> <NA> <NA> <NA> <NA> <NA>
## [61] <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## Levels: 8 9

metadata_Control_V8V9$visit_sum

##  [1] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [10] 8-10    8-10    Control Control Control Control Control Control Control
## [19] Control Control Control Control Control 8-10    8-10    8-10    8-10   
## [28] 8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10    8-10   
## [37] 8-10    8-10    8-10    8-10    Control Control Control Control Control
## [46] Control Control Control Control Control Control Control Control 8-10   
## [55] Control Control Control Control Control Control Control Control Control
## [64] Control Control Control Control Control Control Control Control Control
## [73] Control Control
## Levels: 8-10 Control

metada <- metadata_Control_V8V9%>%
  mutate(visitControl_V8V9 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V8V9 =metada$visitControl_V8V9, subset(metada,select = -c(visitControl_V8V9))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V8V9, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V8V9))
taxtable <- as.data.frame(tax_table(psControl_V8V9))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Phylum)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable,select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)#randomVar = list("+ (1|id)", c("id")), 

raw_p_Control_V8V9 <- metad[1]
corr_p_Control_V8V9 <- metad[2]
effect_size_Control_V8V9 <- metad[3]
status_Control_V8V9 <- metad[4]

psControl_V1 <-
 phyloseq::subset_samples(ps_stool_relab_phylum , visit_cal_cor  %in% c(NA, "1"))

metadata_Control_Control <- as(sample_data(psControl_V1),"data.frame")

metada <- metadata_Control_Control%>%
  mutate(visitControl_V1 = (as.numeric(visit_sum))-1)

metada<- cbind(visitControl_V1 =metada$visitControl_V1, subset(metada, select = -c(visitControl_V1))) # brings visit into the first column

metada <- metada%>%
  dplyr::select(visitControl_V1, id)

# run metadeconfoundR
features <- as.data.frame(otu_table(psControl_V1))
taxtable <- as.data.frame(tax_table(psControl_V1))

#create two-column-dataframe containing corresponding "human-readable" names to the "machine-readable" feature names used as row.names in metaDeconfOutput.  
taxtable <- taxtable%>%
  dplyr::select(Phylum)
taxtable$ASV <- row.names(taxtable)
taxtable<- cbind(ASV=taxtable$ASV,subset(taxtable, select = -c(ASV)))

### run it
metad <- MetaDeconfound(featureMat = features, metaMat = metada,  nnodes=9)

raw_p_Control_V1 <- metad[1]
corr_p_Control_V1 <- metad[2]
effect_size_Control_V1 <- metad[3]
status_Control_V1 <- metad[4]

raw_p <- purrr::map2(raw_p_ControlV2,raw_pControl_V3V5,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V6V7,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V8V9,cbind)
raw_p <- purrr::map2(raw_p,raw_p_Control_V1,cbind)

raw_p_df <- bind_rows(raw_p)
raw_p_df <- data.frame(raw_p_df$Ps)
raw_p_df  <- raw_p_df %>%
  rownames_to_column()%>%
  mutate(p_Control_V2=visit_ControlV2)%>%
  mutate(p_Control_V3V5=visit_Control_V3V5)%>%
  mutate(p_Control_V6V7=visit_Control_V6V7)%>%
  mutate(p_Control_V8V9=visitControl_V8V9)%>%
  mutate(p_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

corr_p <- purrr::map2(corr_p_ControlV2,corr_pControl_V3V5,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V6V7,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V8V9,cbind)
corr_p <- purrr::map2(corr_p,corr_p_Control_V1,cbind)

corr_p_df <- bind_rows(corr_p)
corr_p_df <- data.frame(corr_p_df$Qs)
corr_p_df  <- corr_p_df %>%
  rownames_to_column()%>%
  mutate(q_Control_V2=visit_ControlV2)%>%
  mutate(q_Control_V3V5=visit_Control_V3V5)%>%
  mutate(q_Control_V6V7=visit_Control_V6V7)%>%
  mutate(q_Control_V8V9=visitControl_V8V9)%>%
  mutate(q_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_size <- purrr::map2(effect_size_ControlV2,effect_sizeControl_V3V5,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V6V7,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V8V9,cbind)
effect_size <- purrr::map2(effect_size,effect_size_Control_V1,cbind)

effect_size_df <- bind_rows(effect_size)
effect_size_df <- data.frame(effect_size_df$Ds)
effect_size_df  <- effect_size_df %>%
  rownames_to_column()%>%
  mutate(d_Control_V2=visit_ControlV2)%>%
  mutate(d_Control_V3V5=visit_Control_V3V5)%>%
  mutate(d_Control_V6V7=visit_Control_V6V7)%>%
  mutate(d_Control_V8V9=visitControl_V8V9)%>%
  mutate(d_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

status <- purrr::map2(status_ControlV2,statusControl_V3V5,cbind)
status <- purrr::map2(status,status_Control_V6V7,cbind)
status <- purrr::map2(status,status_Control_V8V9,cbind)
status <- purrr::map2(status,status_Control_V1,cbind)

status_df <- bind_rows(status)
status_df <- data.frame(status_df$status)
status_df  <- status_df %>%
  rownames_to_column()%>%
  mutate(status_Control_V2=visit_ControlV2)%>%
  mutate(status_Control_V3V5=visit_Control_V3V5)%>%
  mutate(status_Control_V6V7=visit_Control_V6V7)%>%
  mutate(status_Control_V8V9=visitControl_V8V9)%>%
  mutate(status_Control_V1=visitControl_V1)%>%
  mutate(ASV=rowname)#%>%
  #dplyr::select(-c(1:10))

effect_table_phylum <- raw_p_df%>%
  full_join(corr_p_df, by="ASV")%>%
  full_join(effect_size_df, by="ASV")%>%
  full_join(status_df, by="ASV")%>%
  full_join(taxtable, by="ASV")

# dplyr::select the entries which have OK_nc in status
effect_table_sig_phylum <- effect_table_phylum%>%
  filter(status_Control_V1=="OK_nc"|status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc")

#pivot long format

effect_table_sig_long_phylum <- effect_table_sig_phylum%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p_"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Phylum, "(Phylum)"))%>%
  dplyr::select(-Phylum)

effect_table_sig_long_phylum$comparison_p <- factor(effect_table_sig_long_phylum$comparison_p, levels = c("Control_V1","Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_phylum%>%
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = as.factor (sign (effectSize)), size = abs (effectSize), color=(fdr))) +
  scale_shape_manual (values = c (25, 24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

# dplyr::select the entries which have OK_nc in status with all differences to Control

effect_table_sig_control <- effect_table_phylum%>%
  filter(status_Control_V2=="OK_nc"|status_Control_V3V5=="OK_nc"|status_Control_V6V7 =="OK_nc" | status_Control_V8V9=="OK_nc"| status_Control_V1=="OK_nc")

#pivot long format

effect_table_sig_long_control_phylum<- effect_table_sig_control%>%
  pivot_longer(cols = starts_with("status"), names_to = "comparison_status", values_to = "status")%>%
  separate(comparison_status, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_status=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  pivot_longer(cols = starts_with("p_"), names_to = "comparison_p", values_to = "raw_p")%>%
  separate(comparison_p, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_p=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_status)%>%
  pivot_longer(cols = starts_with("d"), names_to = "comparison_effectSize", values_to = "effectSize")%>%
  separate(comparison_effectSize, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_effectSize=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_effectSize)%>%
  pivot_longer(cols = starts_with("q"), names_to = "comparison_q", values_to = "corr_p")%>%
  separate(comparison_q, c("variable","Visit1", "Visit_Sum"), "_")%>%
  mutate(comparison_q=paste(Visit1,Visit_Sum, sep="_"))%>%
  dplyr::select(-c(variable, Visit1, Visit_Sum))%>%
  filter(comparison_p==comparison_q)%>%
  dplyr::select(-c(comparison_q, comparison_effectSize,comparison_status))%>%
  mutate(fdr= as_factor(case_when(corr_p <= 0.05 ~ "*", corr_p <= 0.01 ~ "**", corr_p <= 0.001 ~ "***", corr_p <= 0.1 ~ ".")))%>%
  mutate(taxa= paste(Phylum, "(Phylum)"))%>%
  dplyr::select(-Phylum)

effect_table_sig_long_control_phylum$comparison_p <- factor(effect_table_sig_long_control_phylum$comparison_p, levels = c("Control_V1", "Control_V2", "Control_V3V5","Control_V6V7","Control_V8V9"))

effect_table_sig_long_control_phylum%>%
  mutate(shape_sign= as.factor(sign(effectSize))) %>% 
  filter(shape_sign!=0) %>% #fiter out effect sizes that have no direction
  ggplot(aes (x = comparison_p, y = taxa))+ 
  geom_point (aes (fill = effectSize, shape = shape_sign, size = abs (effectSize), color=(fdr)))+
  scale_shape_manual (values = c (25,24)) + 
  scale_fill_gradient2 (low = "blue", high = "red", mid = "white", midpoint = 0) +
  scale_color_manual(values=c("gray22","gray1","gray85"))+
  geom_text (aes (label = stars.pval (corr_p)))+
  theme_grey() +
  theme(axis.text.x = element_text(), #angle = 0, hjust = 0, vjust = 1
       legend.position = "right",axis.text.y = element_text(face = "italic"))+
  scale_x_discrete(labels=c("1", "3", "6-12", "15-18", "21-24"))+
  labs(x= "months from treatment start")

Combine different tax levels into one plot

# combine output tabels
# Select columns from each dataframe
effect_table_sig_long_control <- effect_table_sig_long_control %>%
  dplyr::select(ASV, status, raw_p, comparison_p, effectSize, corr_p, fdr, taxa)

effect_table_sig_long_control_family <- effect_table_sig_long_control_family %>%
  dplyr::select(ASV, status, raw_p, comparison_p, effectSize, corr_p, fdr, taxa)

effect_table_sig_long_control_class <- effect_table_sig_long_control_class %>%
  dplyr::select(ASV, status, raw_p, comparison_p, effectSize, corr_p, fdr, taxa)

effect_table_sig_long_control_phylum <- effect_table_sig_long_control_phylum %>%
 dplyr::select(ASV, status, raw_p, comparison_p, effectSize, corr_p, fdr, taxa)

# bind them
effect_table_sig_all_control <- rbind(effect_table_sig_long_control, effect_table_sig_long_control_family, effect_table_sig_long_control_class, effect_table_sig_long_control_phylum)

# no significant changes on order level observed;  effect_table_sig_long_control_order gets not included here

# retrieve table for supplement
effect_table_sig_all_control_write <- effect_table_sig_all_control %>% 
  mutate(taxa=gsub("ASV\\d{3}_", "", taxa))

write_csv(effect_table_sig_all_control_write, "/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/data/effect_size_table_stool_ControlTimepoints.csv")

#writexl::write_xlsx(effect_table_sig_all_control_write, "/Users/rebecca/Documents/Forschung/IMMProveCF/Manuscript_16S/Submission_CHM/SuppMaterial/SuppTable_18.xlsx")

effect_table_sig_all_control %>%
  mutate(taxa=gsub("ASV\\d{3}_", "", taxa)) %>% 
  arrange(ASV) %>%
  ggplot(aes(x = comparison_p, y = taxa)) + 
  geom_point(
    aes(fill = effectSize, shape = as.factor(sign(effectSize)), size = abs(effectSize), color = fdr)
  ) +
  scale_size_continuous(name = "Absolute effect size (Cliff's delta)", guide = guide_legend(ncol = 5)) + 
  scale_shape_manual(name = "", labels = c("decreased","increased"), values = c(25,24),guide = guide_legend(ncol = 2)) + 
  scale_fill_gradient2(
    name = "Effect size (Cliff's delta)",
    low = "#746999", high = "#69b3a2", mid = "white", midpoint = 0) +
  scale_color_manual(
    name = "", labels = c(". fdr < 0.1", "* fdr < 0.05", "non-significant"),
    values = c("gray22", "gray1", "gray85"), guide = guide_legend(ncol = 3)
  ) +
  geom_text(aes(label = fdr)) +
  theme_grey() +
  theme(
    axis.text.x = element_text(size = 17),
    axis.title.x = element_text(size = 17),
    legend.position = "right",
    axis.text.y = element_text(face = "italic", size = 16),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 17),
    legend.key = element_rect(size = 8),
    legend.key.size = unit(1.5, 'lines')
  ) +
  scale_x_discrete(labels = c("1", "3", "6-12", "15-18", "21-24")) +
  labs(x = "months from treatment start", y = "")

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 110 rows containing missing values or values outside the scale range
## (`geom_text()`).

#ggsave("/Users/rebecca/Documents/Forschung/IMMProveCF/R_analysis/figures/relAb_cuneiform_plot_stool_Controls.png", width = 11, height = 5)